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ABSTRACT 


This study examines various computational techniques to analyze dynamic response and 
failure of sandwich composite materials subject to fluid-structure interaction character¬ 
ized by an acoustic field or the propagation of velocity potential according to the wave 
equation. A displacement-only plate finite element is developed and implemented using 
Discontinuous Galerkin (DG) methodology; its accuracy compares favorably to both the¬ 
ory and Continuous Galerkin methods. Several approaches to analyzing debonding failure 
between skin and core layers of sandwich composite structures are demonstrated and eval¬ 
uated; partial disconnection between neighboring elements at a debonding site shows good 
qualitative agreement with known physical phenomena. A hybrid Finite Element-Cellular 
Automata (FE+CA) approach to modeling an acoustic field with non-reflecting boundary 
conditions is presented, validated numerically, and favorably compared with experimen¬ 
tal results. The FE+CA fluid model is then combined with the DG structural model to 
simulate fluid-structure interaction; this combined model compared favorably with experi¬ 
mental results for the strain field of laminated plates subject to low-velocity impact. Each 
technique addressed shows promise for flexible and accurate modeling of failure initiation 
and propagation in sandwich and laminate composites subject to fluid-structure interaction 
with moderate computational costs. 
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I. INTRODUCTION 


A. MOTIVATION AND OBJECTIVE 

Composite materials can have very beneficial properties for structural applications. 
In particular, polymer composites generally have low density but high strength and stiff¬ 
ness, or a very high specific strength and stiffness. Their resistance to corrosion is another 
tremendous advantage. Consequently, composite materials are used in a number of both 
civil and military applications. Aerospace structures are among the major applications of 
polymer composite materials, especially carbon fiber composites. Increasingly, composite 
materials are used in marine structures, making their use in Naval applications likewise 
more common. 

One of the major difficulties associated with design and analysis of composite struc¬ 
tures is their anisotropic material properties and complex failure modes and mechanisms 
when compared to metallic structures. Anisotropic composite material behavior can be tai¬ 
lored to optimize their use in structures; however, diverse and inter-connected failure modes 
in composites remain a challenge for the composite community. Traditional metallic struc¬ 
tures differ from composites in one other key aspect: resilience and ductility. That is, their 
ability to deform and recover in the elastic zone gives the designer a safety margin much 
larger than that found in the composite world. Composite materials have long been used 
as primary hull materials for sailboats and other small craft, but the nature of larger Naval 
applications of composite materials demands a fuller understanding of the survivability and 
mission impacts of using such structures in the marine environment. 

A great deal of research in FSI has been conducted for aerospace applications, tend¬ 
ing to focus on the effect of the structure on the fluid field. For marine applications with 
a polymer composite structure in contact with water, the comparable density of the com¬ 
posite materials to that of water results in a significant hydrodynamic mass effect on the 
structure. Therefore, the goal of this study is to develop computational techniques to model 
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and simulate transient dynamic responses and failures of composite structures with FSI. To 
this end, it is necessary to develop computational techniques for composite structures as 
well as the fluid medium. Eventually, a method to couple both structural and fluid solvers 
is developed. 

B. PRIOR WORK 

The Finite Element Method (FEM) is a well established tool for generating numer¬ 
ical solutions to the partial differential equations that describe a wide range of physical 
phenomena. In structural solvers, the method combines the geometry and constraints of 
the structure with the material properties of its components to generate a response (e.g., 
displacement, stress, and/or strain) to given loading. This is accomplished by treating the 
structure as a collection of smaller domains (finite elements) of relatively simpler geome¬ 
try, applying and solving the simpler problem on a smaller scale with constraints dictated 
by internal compatibility, and then combining the smaller solutions into a global solution. 
This collection of potentially tedious but simple calculations is an ideal job for a com¬ 
puter. The method of weighted residual used in most finite element codes is a Galerkin 
method in which the test function is the derivative of the trial function with respect to the 
unknown variable(s). Requiring continuity between neighboring elements makes physical 
sense and is the norm, but Discontinuous Galerkin (DG) methods have become an active 
area of research in the last several decades. 

DG approaches to solving boundary-value problems have their genesis in the work 
of Nitsche [1] who first proposed weak enforcement of boundary conditions. Douglas 
and Dupont [2], Arnold [3], Baker [4], and Wheeler [5] expanded the concept to weak 
enforcement of continuity between elements and applied it to elliptic problems. Easing 
of the continuity constraints between elements suggests a number of potential advantages 
to this methodology: element-wise computations lend themselves to parallel computing, 
the independence of the elements from each other allows different orders of interpolation 
within neighboring elements, and the potential to model failure without having to re-mesh 
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is intriguing. In recent years discontinuous approaches have been increasingly applied to 
elliptic partial differential equations such as those that dominate linear elasticity. Arnold 
and colleagues have produced two excellent overviews of the development and characteris¬ 
tics of various methods in [6] and [7]. In particular, they gathered the various formulations 
and cast them all in both flux and primal forms to more clearly see the differences and 
identifying characteristics. Castillo [8] conducted a cost and performance analysis of three 
of these methods. Brezzi, Cockbum, Marini, and Suli [9] discuss the stabilization mecha¬ 
nisms necessary for effective DG formulations. Specific work in DG for solid mechanics 
can be found in the List of References, [10] - [22] among many others. 

Not all parts of a finite element analysis are of equal interest. A homogeneous beam 
or plate, for example, can be modeled with relative fidelity with just a few elements; while 
the analysis of air-flow over a golf ball would require elements of a size comparable to the 
ball’s dimples, which could lead to an exorbitant computational cost for a domain on the 
order of three ball diameters. The obvious solution to such a dilemma is to selectively refine 
the mesh in the area of interest while using a mesh as coarse as possible in other parts of the 
domain. This work seeks to explore the potential of using Discontinuous Galerkin (DG) 
finite element methodology to enable failure to initiate in its natural location anywhere in 
the domain. 

In analysis of structures consisting of composite materials, multiple scales of anal¬ 
ysis can be used in conjunction with each other. That is, the analysis can proceed from 
the micro- (atomic or molecular) level to the meso- (material) level and on to the macro- 
structural) level and back down again as necessary with smearing or decomposition of 
properties and loads as dictated by the analysis required. The properties of a laminated 
fibrous composite, for example, are a function of the properties of both fiber and matrix as 
well as the weave pattern of individual laminae and lay-up pattern of the assembled lami¬ 
nate. Similarly, failure of the same laminate can result from separation between laminae, 
separation of fibers from matrix, or failure of individual fibers. 
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To accurately model FSI an appropriate fluid model is also required. In this work no 
attempt is made to model or solve the full Navier-Stokes equations, a subset that will allow 
compatibility between fluid and solid regimes and approximate the hydrodynamic pressure 
or acoustic field will suffice. Olson and Bathe [23] developed a directly coupled formula¬ 
tion that solves for the velocity potential and hydrostatic pressure in the fluid domain and 
displacements in the structural; this will be the starting point for development of the fluid 
model in this work. In order to apply the fluid model to a maritime domain, appropriate 
boundary conditions must also be included. Two general approaches are to model a vast 
domain and concern oneself with a small subset relatively far from simple but inaccurate 
boundaries-a computationally expensive proposition-or to model an appropriately sized 
domain with non-reflecting boundary conditions-a challenging proposition that remains an 
active area of research [24]. Application of Cellular Automata (CA) to modeling the acous¬ 
tic field following from the work of Chopard [25, 26], Krutar et al. [27], and Kwon and 
Hosoglu [28] will be explored. 

Increasing application of composite technology in maritime applications requires 
further work examining the mass effects imparted to composite structures in contact with a 
fluid environment. In particular, evaluation of damage and residual strength are necessary 
for a proper evaluation of the survivability potential of a proposed design. 

C. ORGANIZATION 

The balance of this thesis is organized as follows: Chapter II includes a short re¬ 
view of linear elasticity followed by development of both full three-dimensional and plate 
elements using DG techniques. Chapter III includes further discussion of the advantages 
and challenges of composite materials and sandwich plates, validation of the developed 
DG structural model’s applicability to sandwich constructs, followed with an examination 
of various failure models. Development and validation of the fluid model incorporating 
velocity potential formulation and non-reflecting boundary conditions comprises Chapter 
IV. Chapter V contains a demonstration of the coupling of the fluid and structural models 
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as well as comparisons with experimental work [29]. Conclusions and recommendations 
for continuing research in this vein will be addressed in Chapter VI. 
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II. STRUCTURAL MODEL 


This chapter will discuss the formulation and various components of the structural 
models used in this work. After first reviewing the equation(s) to be solved and traditional 
numerical approaches, Discontinuous Galerkin (DG) methods will be discussed in general 
followed by a detailed formulation for a three-dimensional solid element. That approach 
will be modified to model a plate element. A discussion of the sensitivities of each element 
to discretization and penalty parameters will ensue. 

A. LINEAR ELASTICITY 

In Voigt notation, the governing equation for the static structural model is the equa¬ 
tion of equilibrium, 

V • <7 + f — 0 (1) 

in which 

® { x &y Try 7~y Z T xz | (2) 

is the stress vector and 

/= {/. fy M T ( 3 ) 

is the body force vector. Combining the constitutive equation of generalized Hooke’s Law, 

a = [D]e (4) 

where e = {s x e y e z ^ xy ^ yz ^ XZ } T is the strain vector and [D] is the 6x6 material 
property matrix, and the kinematic strain-displacement relationship, 

£ = ^(W+W r ) (5) 
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in which u = { Ux u y u z } T = {u v w} T • Equation (1) can be solved for dis¬ 
placement or deflection of a body subject to a load described by f(x, y, z ) and appropriate 
boundary conditions. Displacement at specified points or nodes in the domain of interest 
are the unknown quantities or degree(s) of freedom (dof). 

The desired solution, u(x,y,z), is expressed in terms of the displacement vector 
at each nodal point in the domain and the interpolation functions within each element. If 
we define U as the vector of nodal displacements with the three components for each node 
grouped together, 

u{x,y,z) = N{x,y,z)U (6) 

where 

Hi 0 0 H-2 0 0 . H 8 0 0 

N = 0 Hi 0 0 if 2 0 . 0 H 8 0 (7) 

0 0 Hi 0 0 if 2 . 0 0 H 8 

for a three-dimensional hexahedral element where H % is the three-dimensional cardinal 
basis function corresponding to node i. Similarly, the strain matrix [ B] is defined such that 


e(x,y,z) = B(x,y,z)U (8) 


dHi 

dx 

0 

0 

dm 

dx 

0 

0 ... . 

ym 

dx 

0 

0 

0 

dm 

dy 

0 

0 

dm 

dy 

0 ... . 

. 0 

dm 

dy 

0 

0 

0 

dm 

dz 

0 

0 

dm 

dz ■ ■ ■ ■ 

. 0 

0 

dm 

dz 

dHi 

dm 

0 

dm 

dm 

0 ... . 

dm 

dm 

0 

dy 

dx 

dy 

dx 

dy 

dx 

0 

dm 

am 

0 

dm 

dm 

. 0 

dm 

dm 

dz 

dy 

dz 

dy 

dz 

dy 

am 

0 

am 

dm 

0 

dm 

dm 

0 

dm 

dz 

dx 

dz 

dx 

. dz 

dx 


V, z) = [D]B(x , y, z)U. (10) 
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The dimensions of these elemental matrices and vectors for tri-linear three-dimensional 
hexahedral elements are [U] = 24 x 1 (8 nodes, 3 degrees of freedom per node), [N] = 
3 x 24, [B] — 6 x 24, and [D] — 6x6. U can then be post-processed to calculate the stress 
and strain fields as necessary. Figure 1 shows a notional hexahedral finite element of the 
sort used here. 



1 -1 

Figure 1: Canonical hexahedral finite element 



B. DISCONTINUOUS GALERKIN FORMULATION 

This section will detail the formulation of nodal DG solid (three-dimensional) and 
plate elements to be used to solve for the displacement field in a given linearly elastic 
structure. The plate element to be developed is based on Reissner-Mindlin theory but all 
dof are displacements. 

1. Solid Element 

Liu, Wheeler, and Dawson [21] proposed and implemented a nodal DG formulation 
that can be simply coupled with existing codes for continuous models and can be switched 
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between three different specific methods via the selection of a single scalar parameter, d DG - 
Their formulation is the starting point for this work and is summarized below. 

The domain, f) C R 3 , has a boundary, dQ, comprised of non-intersecting Dirichlet, 
T u , and Neumann, T t , boundaries upon which displacement, u e H l (T u ), and surface 
traction, t e L 2 (T t ) are specified, x = {£j, E 2 , ..., E N } is a non-degenerate discretization 
of Q; in this work, E :J are hexahedra. S = S t + T n T T, is the set of faces of x where Si 
are interior faces. 

Following standard weighted residual methodology, Equation (1) is multiplied by a 
test function and integrated over the domain with the intent of finding a solution that leaves 
no residual from that integral. Since this is a Galerkin formulation, the test function to be 
used is the interpolation function used in the discretization of the domain; as a Discon¬ 
tinuous Galerkin formulation the test function need only be defined within each element. 
Symbolically, 

f a(u) : S/vdV — f (an) ■ vdS — f f ■ vdV. (11) 

J E JdE J E 

Taking advantage of the fact that 

a(u) : Vn = a(u) : Vv T = a(u ) : e(v) (12) 


the elemental equation in question is 


/ a(u) : e(v)dV — / (an) ■ vdS —If' vdV. 

IE JdE Je 


(13) 


Or, summed over all elements, 


: £ ( v ) av - ■ vdS = 1, f ■ vdV. 

E&x^ E dE£S ddE E&x^ E 


(14) 


In a continuous formulation the second term of the last equation would disappear on inter- 
elemental boundaries and only exist on the domain boundary; for DG, however, that is not 
the case. Define jump and average functions across such an interior boundary, between two 
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elements arbitrarily labeled as “Left” and “Right” as 


M = W L - W R 

1 

W = 2 ( w l + w r ). 

Combining those definitions with the identitiy 

[M = MM + MM 


(15) 

(16) 


(17) 


and the assumption that traction across interfaces is continuous, results in 



E 

dE&{S-r u )' 


'dE 


{a(u)n s } ■ [ujdS 1 


Y [ f-vdV+ Y [ t-vdS. (18) 

E&x E dEeF't ^ dE 


Liu et al. [21] add face integrals J gE {a(v)n s } ■ [u]dS and Mp f 9E [u] ■ [v}dS, both of 
which disappear for an exact solution, to control symmetry and provide stabilization; they 
generate the following bilinear and linear forms: 



a(u) : e(v)dV — 'Y] 

dE£{Si+ r„) 


+ 'Y, ® dg 

dE&(Si+T u ) 



{a(u)n s } ■ [ujdS 1 


JdE 

{a(v)n s } ■ [wjdS' 


+ £ 

dE&(Si+T u ) 


8 P G 

Isl 


IdE 


u\ • [u]d£ (19) 
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L ( v ) = [ f ■ vdV + ^ f t ■ vdS + ^ 9 DG f {cr(v)n s } • udS 

E&x E 9E£ r t ^ dE OEgFu J dE 

+ Y'' I u ■ [t^dS 1 (20) 

dEe r u I s ! ^ dE 

where G is the shear modulus, 5 P is a scalar penalty parameter, \s\ is the square root of 
the area of the element’s face, and 6 DG indicates the DG method in use. For non-zero 5 P , 
if 9og = — 1 the method is the Symmetric Interior Penalty Galerkin (SIPG), which also 
corresponds to the Local Discontinuous Galerkin (LDG) method of Cockburn and Shu with 
(3 = 0 [30]; if 9pc = +1 it is Non-symmetric Interior Penalty Galerkin (NIPG); if 9 DG = 0 
it is Incomplete Interior Penalty Galerkin (IIPG). If 9dg = 0 and d p = 0 the method is that 
of Oden, Babuska, and Baumann (OBB). The problem statement is now: find u G V such 
that 

a(u,v) = L(v) VveV (21) 

H\x) = {v e L 2 (Q) : v\ Ej G H^E^Ej Gy}; V — {v G H\ X )}- (22) 

Equations (19) and (20) can be converted into matrix-vector form. The integrand of 
the first term of (19) is: 

a(u ) : e(v) = DB(x,y, z)U : B(x,y,z) = B 7 (x,y, z)DB(x,y, z)U. (23) 

Since U is independent of (x, y, z), it can be taken outside the integral 

[ a(u) : e(v)dV = [ B T DBdVU = K V U (24) 

J E J E 

so that the matrix, K v resulting from the volume integral is multiplied by the displacement 
vector. This K v is the identical to the elemental stiffness matrix common to continuous 
Galerkin formulations. 

In a similar fashion, the surface integrals of Equation (19) can also be expressed 
as matrix-vector products. Once expanded the vector components will be comprised of 
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various products of the unknown vector, u L or u R and the test function, v L or v R . The 
interface stiffness matrices will be referred to as K l MN where the superscript, i refers to the 
integrals in the order they appear in Equation (19); the subscripts will take the values L and 
R, referring to the left and right sides of the interface, respectively. The subscript M will 
correspond to the v component and the N will correspond to the u component. 

For example, the expansion of the first surface integral into the K\ IN components 
is 


{cr(u)n s } ■ [v]dS = - f + a ^ UR \ s ) . (y L - VR )dS 

> Js z 

= f {v{uL)n s ) ■ v L dS + ^ f ( a{u L )n s ) ■ v R dS 
J S J s 

- \ f ( a(u R )n a ) • v L dS + \ f ( cr(u R )n s ) • v R dS. (25) 


Expressing the normal vector in terms of a matrix of the direction cosines, 


A s = 


n x 0 0 ri y 0 n z 


0 n y 0 n x n z 0 
0 0 n z 0 n y n x 


(26) 


the terms of Equation (25) can be expressed as 


2 [Wu L )n s ) • v L dS =~\f N T L A s DB L dSU L = K\ l U l 
S J s 


(27) 


2 J{<?{u L )n s ) • v R dS = 2 ^, N T R A s DB L dSU L = K X RL U L (28) 

- ^ f{<u R )n s ) ■ v L dS = ~\J s N T L A s DB R dSU R = K\ r U r (29) 

J [ (v(u R )n s ) ■ v R dS = \ [ N T R A s DB R dSU R = K l RR U R . (30) 
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Conveniently the interface stiffness matrices that result from the third term of (19) are 
simple re-arrangements of those calculated from the second term. That is, K\ h = Kjj', 
K RR = I\ RR r , K\ r = K rj T , and K R[ = K\ R . The final term of (19) is referred to as 
the interface penalty stiffness, K\j N and its components are calculated as follows: 


SpG 

Isl 



u] ■ [u]dS' 


S P G 


(u L - u R )[y L - v R )dS 


( u L v L - u L v R - u R v L + u R v R )dS (31) 


1*1 Js 

Yr [ u L v L dS = [ NlN L dSU L = K\ l U l (32) 

\s\ Js |s| Js 

- 5j fr I u L v R dS = [ N T R N L dSU L = K\ l U l (33) 

l s l Js l s l Js 

- Yr f u ^ dS = ~TT f N T L N R dSU R = K 3 lr U r (34) 

M Js ! s l Js 

^ [ u R v R dS = f N T R N R dSU R = K rr U r . (35) 

l s l Js l s l Js 


The terms of L(v) are vectors resulting from integrating the body forces and given 
boundary conditions (both displacement and traction) over the volume and surfaces. 


/ / • vdV — / N 1 fdV = F b 

' E J E 


where / = (f x , f y , f z ) T for the element in question. 


/ t ■ vdS — / N T tdS = F t 

'dE JdE 


/ {a{v)n s } ■ udS = / B 1 D 1 (A 8 ) 1 udS = F' u 

'dE JdE 


Yv [ [v]dS = Y f 
l s l JdE l s l JdE 

\T „.,\T 


N T udS = F p 


where t = (t x , t y , t z ) T and u = (u, v, w) T for the boundary face in question. 


(36) 


(37) 

(38) 

(39) 
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To ensure each interface integral is calculated once and only once, as the loop 
through the elements progresses surface integrals are calculated only for those faces corre¬ 
sponding to the positive canonical directions. This convention relies on element numbering 
also proceeding in the positive directions, or element e’s +r neighbor is an element num¬ 
bered greater than e. 

One of the advantages of discontinuous Galerkin formulations is the independence 
of the elements with the exception of the numerical fluxes between them. If this were 
a time-dependent problem or if there were some other source of data about neighboring 
elements, the various terms of a(u, v ) could be computed for each element using the dis¬ 
placement data from the previous time step as the Ur terms for its neighbors. In a static 
treatment, the various K l MN matrices will be assembled into a total system stiffness matrix; 
likewise, the various F l vectors will be assembled into a total system force vector and a 
global KU = F and solved simultaneously. 

Assembly of the global [K] matrix highlights the most obvious difference between 
CG and DG methods: the relationships between the elements. Figures 2 and 3 illustrate 
the relationships between neighboring elements in the two methods. The interface between 
the DG elements may have the same geometric location initially, but the faces are treated 
as separate entities; each element’s dofs are wholly contained in that element. The CG 
elements actually share the face/edge and have common nodes and dofs there. Figures 4 
and 5 show the connectivity and sparsity patterns for the assembled global stiffness matrices 
of the same four element by two element by one element discretization generated by the 
two methods. Cursory examination of Figure 5 reveals the elemental connectivity: the full 
blocks along the diagonal are the eight elemental stiffness matrices; the sparser off-diagonal 
blocks represent the interfaces. It is apparent, therefore, that element one’s neighbors are 
element two and element five. The connectivity of Figure 4 is related to the individual 
nodes and dofs as opposed to the elemental connectivity of Figure 5. These two figures 
also illustrate one of the significant disadvantages of DG methods - a tremendous growth 
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in the scope of the numerical problem; these two matrices can be used to solve the same 
physical problem, yet the CG version is only half as large. 


Q- 




x 


0 - 




O left 
□ common 
^ right 


x 


Figure 2: Connectivity between adjacent CG elements. After [31] 
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Figure 3: Connectivity between adjacent DG elements. After [31] 


CG 4x2x1 



0 20 40 60 80 


nz = 2979 


Figure 4: Sparsity pattern of CG global stiffness matrix for a 4x2x1 element structure 
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DG 4x2x1 



0 50 100 150 

nz = 8128 


Figure 5: Sparsity pattern of DG global stiffness matrix for a 4x2x1 element structure 
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a. Validation 

A pair of beam problems were used to validate code derived from this 
formulation: a cantilever subjected to a concentrated force at the free end and a simply- 
supported beam subjected to a concentrated force at mid-span. Both beams are 8m long 
and have unit cross-sectional area and the same isotropic material properties. Errors relative 
to maximum deflections predicted by Euler Beam theory were calculated for various dis¬ 
cretizations and plotted in Figures 6-7, both of which demonstrate quadratic convergence 
rates (the predicted rate for linear elements) for all three methods (NIPG, SIPG, IIPG). 



Figure 6: Convergence of three-dimensional DG formulation vs. discretization for can¬ 
tilever beam deflection 


2. Plate Element 

Adapting the above formulation to different types of elements is a matter of using 
appropriate [N], [B], and 77 matrices. Kwon and Bang [32] developed a displacement 
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Figure 7: Convergence of three-dimensional DG formulation vs. discretization for simply- 
supported beam deflection 


only plate model following Mindlin-Reissner plate theory with the following definitions. 

i I {ttl Vi Ui+ n Vi+ n W\ U2 V2 ri2+n ^'2 • n VJ2 ■ ■ ■ U n V n U 2 n V2n ^n}’ 

(40) 

where n is the number of nodes on the bottom of the plate, node i + n is taken to be above 
node i and transverse deflection of top and bottom are taken to be equal, or iv, = w i+n , 
eliminating the need for w i>n . 


{ e &} - {e x e y 7 xy } T - [Bb\{U} 

(41) 

{(?b} = {<j x cr y r xy } T = [D b ]{e fe } 

(42) 

{e.}= {lyz lxz} T =[Bs]{U} 

(43) 

{<**}= {r yz t xz } t = [77 s ]{eJ 

(44) 
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where the subscripts b and s indicate bending and shear, respectively. The volume integral 
that is used as the stiffness matrix then takes the form 

[K\= [ [B b f[D b ][B b ]dn+ [ [B a ] T [D 8 ][B a ]dn (45) 


[B 

d= [[Bn] [B b2 ] 

[Bb 3 ] 

[Bb 4]] 


(46) 


1 

H | s < 

0 

H2 d ~§ 

0 

0 


[Bu] = 

0 


0 

h 2 ^± 

z dy 

0 

(47) 


si® 

Hi 9 ~§ 

U dNi 
n 2~dy 

U dNi 

0 



[Bs]= [[S sl ] [B s2 ] [B s3 \ [B s 4 ]] (48) 


[Bsi\ = 


0 


0 fi 

OX 


0 


n n dNi 

ill A 

L oz ox 


In the bending and shear strain matrices, [B b ] and | B s \, II, (x. y) are the two dimensional 
nodal interpolation functions in the planar directions and Ni(z) are the one dimensional 
transverse nodal interpolation functions. The shear term must be under-integrated numeri¬ 
cally to prevent shear locking for thin elements. Further discussion of the numerical inte¬ 
gration schemes used can be found in Appendix B. An approach used in other DG formu¬ 
lations of the plate bending problem is to use a lower order interpolant for shear terms than 
that used for bending terms [15]. 


a. Validation 

A clamped square plate subject to a concentrated load applied at its center 
was used to validate the plate element model. The plate is .3048 m x .3048 m x .00635 
m (12 in x 12 in x 1/4 in) and isotropic. Theoretical values were calculated according 
to Timoshenko [33]. All three methods demonstrate quadratic convergence when linear 
elements are used, as shown in Figure 8. 
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Figure 8: Convergence of center deflection of a clamped plate comprised of plate elements 

A second validation was conducted by comparing dynamic continuous and 
discontinuous models of the same square plate. Both models consist of 144 square elements 
such that the characteristic length of each element is four times its thickness. Clamped 
boundary conditions and zero deflection and velocity initial conditions were applied and 
both models have lumped (diagonal) mass matrices. Figures 9 and 10 show excellent agree¬ 
ment between the two models for the transverse displacement and velocity of the center of 
the plate subjected to a constant concentrated force applied at its center. Appendix A details 
the time integration algorithms used. 

One further validation compares the present formulation with central deflec¬ 
tion of uniformly loaded composite plates as described by Lok and Cheng [34]. A square 
discretization with progressively more elements per side in plate and a single thickness 
element were modeled. 
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Figure 9: Comparison of displacement calculated for Continuous Galerkin and Discontin¬ 
uous Galerkin models of a center loaded clamped plate 



Figure 10: Comparison of velocity calculated for Continuous Galerkin and Discontinuous 
Galerkin models of a center loaded clamped plate 


23 










































x 10“ 3 



Figure 11: Convergence of central deflection of uniformly loaded clamped laminated plate, 
comparison with [34] 


24 







3. Stabilization and Penalty Parameters 

One challenge in using these discontinuous Galerkin elements is selection of an 
appropriate penalty parameter for the last two terms of Equations (19) and (20). As 5 P 
approaches infinity the methods return to their continuous roots, so selection of too large 
a penalty is a waste of computing resources. The penalty must also be large enough to 
guarantee existence and uniqueness of the solution [35]. Additionally, the nature of the 
penalty is often described as a function of the local element order or size with little other 
clarification. The nature of penalizing the jump as a stabilization mechanism is discussed 
by Brezzi et al. [9] for elliptic DG formulations in general. Others have computed lower 
bounds on penalties for various formulations, [36], [37], [38] among others. 

In Liu’s formulation [21], the penalty term is a surface integral multiplied by a 
parameter to be determined and the material shear modulus divided by the square root of the 
area over which the integral is calculated. Another intriguing approach is that of Ainsworth 
and Ran ki n [39] in which they compute a lower bound on the penalty parameter that is a 
function of the method selector ( Odg above) and the maximum eigenvalue of the elemental 
stiffness matrix; that value is also divided by a term analogous to the square root of the area 
of the integral. In view of this approach, the presence of the shear modulus in Liu’s penalty 
term serves as a scaling factor to keep the penalty in the same numerical neighborhood 
as that of the volume integral, K v or elemental stiffness matrix. Anticipating using this 
formulation for a composite material with potentially vastly different shear moduli between 
elements, we will combine these two approaches, replacing S p G with 8 P - ma x(p(K v )). As 
long as 6 P > (1 — d DG ) 2 , a unique solution will exist [39]. 
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III. COMPOSITE STRUCTURES 


Composite materials can provide designers with optimal combinations of strength, 
weight, flexibility, and other physical characteristics for their application. Effective design, 
however, requires a thorough understanding of material behavior in the expected operat¬ 
ing environment-data which can be cumbersome to obtain, making the development of 
effective simulations key. 

The strength of a chain is determined by that of its weakest link, but the utility of 
composite materials is the improvement each constituent element brings to the desired ma¬ 
terial properties of the whole. Like alloys, particulate composites generally demonstrate 
better properties than the homogeneous matrix material by virtue of the reinforcement pro¬ 
vided by specifically chosen additives. 

A. MULTI-SCALE MODEL 

Analysis of structures consisting of composite materials, and of the materials them¬ 
selves, require multiple scales of analysis to be used in conjunction with each other. That 
is, the analysis can proceed from the micro- (atomic or molecular) level to the meso- (ma¬ 
terial) level and on to the macro- (structural) level and back down again as necessary with 
smearing or decomposition of properties and loads as dictated by the analysis required. The 
properties of a laminated fibrous composite, for example, are a function of the properties of 
both fiber and matrix as well as the weave pattern of individual laminae and lay-up pattern 
of the assembled laminate. Similarly, failure of the same laminate can result from sepa¬ 
ration between laminae, separation of fibers from matrix, or failure of individual fibers. 
An illustration of this process is shown in Figure 12. Kwon [40] presents an extensive 
discussion of these cycles. 

This study’s focus on structural responses of existing composites allows the ac¬ 
ceptance of micro-level analysis already conducted in the design and construction of the 
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Stiffness Loop 


r -a 



Stress Loop 


Figure 12: Multi-scale analysis cycle for a fibrous composite. From [40] 


composite materials that will be assembled to create the sandwich plates of interest. That 
is, we will not deal with the properties of the fibers that are woven into E-glass, but will 
instead treat the manufacturer’s given orthotropic properties of a ply as known. Specific 
material properties used are listed in Appendix B. 

B. SANDWICH COMPOSITES 

This work will attempt to model sandwich composites for plate and shell applica¬ 
tions. These materials are comprised of low density cores that are relatively stiff trans¬ 
versely and skin layers that provide in-plane strength to the structure. Because the core 
material is used to provide transverse stiffness, it is tempting to model it using full three- 
dimensional solid finite elements; however, the aspect ratio of the plate structure and its 
included elements makes solutions generated using a plate element for all layers more ac¬ 
curate. That is, while the core is generally the thickest component in the sandwich, it is 
still thin relative to its planar cross-section, giving it a sub-optimal aspect ratio for solution 
with three-dimensional elements. This model of a composite material will consist of an as- 
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semblage of elements that are each homogeneously comprised of the constitutent materials 
and of the type of element described. That is, the elemental models developed in Chapter 
II will be used with [D] matrices describing the material properties of constituent materials 
inserted appropriately. Care must be used to ensure the material property matrix contains 
not only the correct individual properties for each element, but also that it is of the cor¬ 
rect form. Modeling an orthotropic material with a [D] matrix calculated in isotropic form 
may result in significant loss of accuracy. Figure 13 illustrates the two different transverse 
lay-ups used to model sandwich plates in this work; the three layer model includes only 
the core and two skin faces; the five layer or “with resin” model includes a relatively thin 
layer with material properties similar to common adhesives used in assembling sandwich 
composites. The two models vary slightly in thickness as well as in overall stiffness due to 
the inclusion of the extra layers. 


Three layer plate model 


skin 


core 


skin 


Five layer plate model 



Figure 13: Three and Five Layer Sandwich Plate models 

Schmit and Monforton [41] developed a discrete element method to predict the 
static deflection of sandwich plates and shells with laminated faces under a variety of 
boundary conditions. Kanematsu and Hirano [42] expanded on that work to examine both 
bending and vibration of sandwich plates; their work also included experimental validation. 
The deflection of a 50 inch square clamped plate with a one inch thick core of aluminum 
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honeycomb faced by two 0.015in thick aluminum skins under uniform pressure was calcu¬ 
lated by both papers and is used here to validate the DG structural model developed in the 
previous chapter. Using one thickness element for each face and one for the core, maxi¬ 
mum static deflection was calculated for a range of square planar discretizations; they are 
plotted in Figure 14. Good agreement was reached with as few as four elements in each 
direction, and refinement further than twelve elements in each direction was shown to be 
unnecessary. Both sets of authors neglected in-plane bending of the core material, which 
the present formulation does not, therefore the current model is sitffer and returns a slightly 
lesser static deflection. 

Relative convergence of both maximum deflection and maximum bending stress 
in a clamped, square, five-layer sandwich plate subject to a concentrated center force was 
examined by modeling a quarter of the plate, taking advantage of symmetry to achieve 
finer discretizations without incurring excessive computational cost. For this study, a range 
of twelve to thirty elements per side of the quarter-plate was modeled. Deflection is the 
primary variable and bending stresses are post-processed quantities. Both skin and core 
stresses were calculated, but as there was minimal difference between the two, core stresses 
are omitted from the plot for clarity. Convergence was calculated relative to the finest 
discretization calculated, thirty elements per quarter-plate side (dx=0.0075 m) and is shown 
in Figure 15. Both quantities converge at better than quadratic rates, the predicted rate for 
linear elements. 
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Figure 14: Convergence of DG Sandwich Plate max deflection (due to uniform load) to 
predicted, comparison with [41, 42] 



Figure 15: Relative convergence of DG Sandwich Plate max deflection and skin stress (due 
to concentrated load) 
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C. FAILURE MODES AND IDENTIFICATION 

Failure of sandwich structures is a function of the constitutent materials, the ge¬ 
ometry of the structure, and the nature of the loading. Common failure modes of such 
structures include debonding, delamination, core crushing, skin wrinkling, and general 
buckling. Debonding is the separation of the skin material from the core and delamination 
usually refers to the separation of layers within the skin material. In this work, debonding 
will be the primary failure mode examined. No attempt is made here to develop failure 
criteria; we will instead attempt to develop an appropriate method to reflect failure within 
a structural model. 

Previous work [43] concluded that modeling an independent layer representing the 
adhesive between laminae of composites was necessary to observe the delamination failure 
mode. In this work we will examine whether such a layer can be omitted when Discontin¬ 
uous Galerkin (DG) techniques are applied in assembling the structure. 

One-quarter of a twenty-four by twenty-four element square clamped plate with the 
same aluminum skins and honeycomb core as was used in the previous section was the 
basis for this examination. This plate model is 450 mm x 450 mm, the core is 10 mm 
thick, each skin is 0.375 mm thick, and the load is a concentrated force of 1000 N applied 
to the center of the plate. The model was assembled as described in Chapter II and the 
global displacement vector was calculated. The degrees of freedom of interest, those on 
the interface between the bottom of the core and the top of the lower skin, were extracted 
and post-processed to calculate the bending stress vector at each point on that interface. 
In order to display the calculated data, an interpolation function describing resulting stress 
values of interest, cr x , was generated using MATLAB’s TriScatteredlnterp function 
and then applied over a grid of the same dimension as the original discretization and plotted 
using contourf. 

Figures 16 and 17 show the resulting normal stress in the x direction on the top 
of the lower skin (the side facing the core) and on the bottom of the core (the side facing 
the lower skin) for a three layer model. The upper right-hand comer of these quarter-plate 
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plots corresponds to the center of the whole plate and is the point of application of the 
concentrated load. As expected, the peak stress values for both components are located at 
the center of the plate, and the stiffer skin is taking a significantly larger portion of this in¬ 
plane load. The effects of the clamped boundary conditions can also be seen at the edges of 
both graphs. This process was repeated for a five-layer model with a lay-up of: skin-resin- 
core-resin-skin. In this case the resultant stresses were calculated for the top of the lower 
skin and bottom of the core as before as well as both the top and bottom of the intervening 
resin layer. The results are displayed in Figures 18 - 21. The stresses on the two faces 
of the resin layer are close enough in magnitude to treat them as equal. Additionally, the 
stresses on the skin and core are negligibly affected by the insertion of the resin layer into 
the model. Unfortunately, the relative magnitudes of the stresses do not suggest an intuitive 
or convenient stress-based failure criteria that could be applied to the resin in absentia 
based on the calculated values on the faces of the skin and core. Therefore, the inclusion 
of an interface layer is needed to accurately model debonding of skin faces from cores of 
sandwich composites. 



Figure 16: Planar stress on the skin of a clamped three-layer sandwich plate subject to 
concentrated center force 
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Figure 17: Planar stress on the core of a clamped three-layer sandwich plate subject to 
concentrated center force 



Figure 18: Planar stress on the skin of a clamped five-layer sandwich plate subject to 
concentrated center force 
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Figure 19: Planar stress on the core of a clamped five-layer sandwich plate subject to 
concentrated center force 



Figure 20: Planar stress on the bottom of a resin layer of a clamped five-layer sandwich 
plate subject to concentrated center force 
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Figure 21: Planar stress on the top of a resin layer of a clamped five-layer sandwich plate 
subject to concentrated center force 
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D. FAILURE MODELING 


This section will compare and contrast three proposed methods of modeling debond¬ 
ing failure between the core and the resin layer opposite an imposed concentrated force at 
the center of a sandwich plate. Each of the methods to follow will be demonstrated using a 
five layer plate discretized into a twenty-four by twenty-four element mesh. The calculated 
planar stress values on the bottom of the core, the top and bottom of the lower resin layer 
and the top of the skin layers will be displayed and discussed. The undamaged model as 
displayed in Figures 18-21 will be used as a baseline reference case. 

1. Damage via Complete Disconnection 

Mergheim et al. [22] introduce a scheme that combines DG methods with existing 
interface methods for modeling failure. Specifically, they propose re-defining Equation 
(19) as follows: 

a(u,v ) = f a ( u ) '■ £(v)dV — ^2 (1 — a) f {cr(w)n s } • [njdS 1 

egx ^ e dE&(Si+ r„) ^ dE 

+ ^2 0dg(1 ~ «) f {cr(v)n s } • [t^dS 

dE&(Si+ r„) JdE 

+ ^2 [ M ' f(l — a )“TT"[ M ] + oT['u]'') dS 1 . (50) 

dEe(Si+r u ) ^ dE ^ I s ! ' 

where a is a switching factor and t is a traction vector governed by a traction-separation 
law. In the pre-critical or undamaged regime, a = 0 and Equation (50) is identical to 
Equation (19); in the post-critical or damaged regime, a = 1 and the surface integrals 
representing connected interfaces are replaced by a traction-separation law that models 
progressive failure. 

Adapting this concept and assuming complete failure of the interface between the 
core and resin, debonding damage was simulated by separating a four by four array of 
resin elements surrounding the center of the plate from their core element neighbors. In 
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the present formulation, all of the area integrals of Equation (19) are calculated for both 
interior and Dirichlet exterior boundaries. If the interface between two elements is deemed 
to have failed, those faces can then be considered members of T u rather than Si, so the 
K)j and K l RR terms remain and the various left and right components {K) R and K' RL ) are 
simply deleted from the global K matrix. 

All other parameters were unchanged from the undamaged case. The resulting 
stress fields are displayed in Figures 22-25. The skin stresses are largely unchanged from 
the undamaged state. The core, however, bears the brunt of this simulated debonding; its 
maximum stress value is twenty times that of its undamaged version. This is because this 
simulation has effectively removed all constraints on/supports to the core in the area of 
greatest load. In the undamaged model, the static deflection at the center of the plate is 
consistent through the thickness; that is, the transverse displacement at all center nodes, 
in both skin layers, both resin layers, and the core, has been the same. In this example, 
the static deflection is consistent from the top of the structure (point of load application) 
down to the bottom of the core; the deflection of the lower resin and lower skin layers was 
consistent within those two layers, but markedly less than that above. Deflection curves 
along the centerline of the plate for the bottoms of the core, resin, and lower skin layers are 
shown in Figure 26. This is because this method of debonding the plate has also eliminated 
any means of transferring the load between those layers within the damage zone. Complete 
disconnection of inter-elemental interfaces is not the proper way to model this sort of dam¬ 
age. A more consistent approach may be to disconnect the planar dofs between elements, 
but leave the transverse dofs connected. 
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Figure 22: Planar stress on the skin of a damaged via complete disconnection clamped 
five-layer sandwich plate subject to concentrated center force 
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Figure 23: Planar stress on the core of a damaged (via complete disconnection) clamped 
five-layer sandwich plate subject to concentrated center force 
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Figure 24: Planar stress on the bottom of the resin layer of a damaged (via complete dis 
connection) clamped five-layer sandwich plate subject to concentrated center force 



Figure 25: Planar stress on the top of the resin layer of a damaged (via complete discon 
nection) clamped five-layer sandwich plate subject to concentrated center force 
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Figure 26: Static deflection along centerline of damaged (via complete disconnection) 
clamped five-layer sandwich plate subject to concentrated center force 
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2. Damage via Partial Disconnection 

To implement a more physically consistent disconnection between core and resin 
elements in the debonding zone, only those rows and columns of the interface sub-matrices 
that correspond to the planar (rq and v,) dofs will be removed. That is, the entries in those 
sub-matrices that correspond to the transverse (wy) dofs will be left in place. Once again, 
this removal of entries from the global stiffness matrix is executed after its assembly - a 
step that can be repeated as necessary for progressive failure with relative simplicity. The 
resultant stresses can be seen in Figures 27 - 30. The maximum stress values are consis¬ 
tently lower than the undamaged case in all three materials, but more noteworthy is the 
movement of the location of maximum stress from the center of the domain, directly below 
the imposed load, to the edge of the debonding zone. This behavior appears to be consistent 
with a physical stress concentration on the edge of a discontinuity in a structure and may 
be useful in modeling damage propagation. Figure 31 shows that the static deflection cor¬ 
responding to this failure model is physically consistent: the core does not deflect beyond 
or through the resin and skin layers below it, but all three layers show sharper deflection 
within the debonding zone than in the rest of the domain. 

Figures 32 - 34 show the stress profile generated using this partial disconnection 
method with the quarter-plate model and a discretization of thirty elements per side. The 
stress concentration effect observed in the relatively coarse meshes of Figures 27 - 30 are 
clearly present in the finer model as well. 


42 




Figure 27: Planar stress on the skin of a damaged (by partial disconnection clamped) five 
layer sandwich plate subject to concentrated center force 



Figure 28: Planar stress on the core of a damaged (by partial disconnection clamped) five 
layer sandwich plate subject to concentrated center force 
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Figure 29: Planar stress on the bottom of the resin layer of a damaged (by partial discon¬ 
nection) clamped five-layer sandwich plate subject to concentrated center force 
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Figure 30: Planar stress on the top of the resin layer of a damaged (by partial disconnection) 
clamped five-layer sandwich plate subject to concentrated center force 
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Figure 31: Static deflection along centerline of damaged (via partial disconnection) five 
layer sandwich plate subject to concentrated center force 



Figure 32: Planar stress on the skin of a damaged (by partial disconnection) clamped five 
layer sandwich plate subject to concentrated center force, fine view 
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Figure 33: Planar stress on the core of a damaged (by partial disconnection) clamped five- 
layer sandwich plate subject to concentrated center force, fine view 



Figure 34: Planar stress in the resin layer of a damaged (by partial disconnection) clamped 
five-layer sandwich plate subject to concentrated center force, fine view 
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3. Damage via Reduced Moduli 

The next method of imposing damage in this model was to leave all elements con¬ 
nected as in the undamaged state, but to reduce the effectiveness of the resin elements in 
the damage zone. The same resin elements identified in the previous attempt remained 
connected to their core counterparts, but their shear and elastic moduli were reduced to 1% 
of the values used for the rest of the layer-an arbitrarily chosen reduction. The resultant 
stresses are plotted in Figures 35 - 38. In this model, the stress in the core elements is com¬ 
parable to that of the undamaged plate and the dramatically lower stress values in the center 
of the resin layer is entirely attributable to the lower modulus. This method would seem to 
be similar to Mergheim’s insertion of a traction-separation law to describe the progression 
from damage initiation to complete separation [22]. 



Figure 35: Planar stress on the skin of a damaged (by reduced modulus) clamped five-layer 
sandwich plate subject to concentrated center force 
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Figure 36: Planar stress on the core of a damaged (by reduced modulus) clamped five-layer 
sandwich plate subject to concentrated center force 



Figure 37: Planar stress on the bottom of the resin layer of a damaged (by reduced modulus) 
clamped five-layer sandwich plate subject to concentrated center force 
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Figure 38: Planar stress on the top of the resin layer of a damaged (by reduced modulus) 
clamped five-layer sandwich plate subject to concentrated center force 
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E. SYNTHESIS 


Table 1 contains a summary of the above results; the numerical stress values dis¬ 
played are of use only in their relationship to each other. These were all calculated from 
the twelve by twelve element quarter-plate, which corresponds to the coarsest discretization 
used in the convergence calculations displayed in Figure 15. 

Complete elimination of interface terms from the global stiffness matrix does not 
correctly model physical constraints on the core from remaining layers in the structure be¬ 
low the section deemed to have been delaminated. Disconnection of the planar terms whilst 
retaining transverse interface terms does appear to correctly model expected physical be¬ 
havior. Adjusting the physical properties of the resin or interface layer remains a potential 
tool for faithful modeling of traction-separation laws. 



Undamaged 

Complete 

Disconnect 

Semi- 

Disconnect 

Reduced 

Modulus 


value 

location 

value 

location 

value 

location 

value 

location 

skin 

96.88 

center 

77.89 

center 

41.02 

zone edge 

97.61 

center 

core 

0.280 

center 

4.389 

center 

0.111 

zone edge 

0.284 

center 

resin top 

10.82 

center 

8.639 

center 

4.413 

zone edge 

4.331 

zone edge 

resin bottom 

10.85 

center 

8.724 

center 

4.641 

zone edge 

4.303 

zone edge 


Table 1: Comparison of maximum planar stress values (in MPa) and locations for damaged 
and undamaged clamped sandwich plates subject to a concentrated center force 
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IV. FLUID MODEL 


Following the work of Olson and Bathe [23], Fluid-Structure Interaction (FSI) will 
be analyzed via a velocity potential in the fluid. The transverse plate velocity will be 
matched to the z component of the fluid velocity as a compatibility condition between the 
two domains. The scalar velocity potential follows the wave equation, so a sufficiently ac¬ 
curate and efficient model of that equation with appropriate boundary conditions is required 
for the fluid portion of this work. 

A great many methods are available to model the fluid mechanics necessary for this 
work, but as our primary interest is in the structural side of the fluid-structure interaction, 
achieving adequate accuracy without incurring significant computational cost lead to the 
exploration of CA methods. 

A. VELOCITY POTENTIAL AND WAVE EQUATION THEORY 

The velocity potential, (f), in the fluid domain is defined as 

v = V< j> (51) 

where v is the velocity of the fluid. The wave equation is 

u = c 2 V 2 u (52) 

where c is the acoustic speed of the (fluid) medium. Coupled with appropriate initial con¬ 
ditions, this well-posed initial value problem has been much studied and discussed [44]. In 
one-dimension, the problem 


Utt = C 2 U XX — OO < X < oo 0 < t < oo 


«(z,0) = f(x) 


(53) 
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U t (x, 0) = g(x) 


is satisfied by the D’Alembert solution 

1 1 r x + ct 

t) = -[f(x - ct ) + f(x + ct )] + — / g(Qd£. (54) 

" Jx—ct 

The three-dimensional extensions are 


u tt = c 2 V 2 u (x, y, z) e R 3 


(55) 


u{x,y,z, 0) = <p(x,y,z) 
u t (x,y,z, 0) = ip(x, y, z) 

which is satisfied by 

d 

u(x,y,z,t) = tip + — [ttp\ (56) 

where ip and (p are the averages of their respective initial conditions over the sphere of 
radius ct centered at (x, y, z); specifically, 


f* 7T p2ir 


x,y,z ) = 


4nc 2 t 2 


'o Jo 


ip(x + ct sin 4b cos 9,y + ct sin 0 sin 9, 

z + ct cos 0) (ct) 2 sin 0d#d0 (57) 


and 


7T p2lT 


<p(x, y, z) = ,o 2 / / <p(x + ctsm<pcos9,y + ctsm(psm9, 

4:70C~t J Q J Q 


z + ct cos0)(cf) 2 sin0d6 ) d0. (58) 


The integrals in (57) and (58) are rarely simple to evaluate analytically. 
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B. FINITE ELEMENT MODEL OF THE WAVE EQUATION 


Beginning with the wave equation as applied to the velocity potential, 

2_ 9 , d 2 (j) 

multiply all terms by a test function, w, and integrate over the domain to get 

f wcp dfl — c 2 f w\/ 2 (p dO = 0 

Jn Jn 

integrate the second term by parts and apply Green’s Identity to get 


/ w(J) dfl + c 2 / Vw ■ VhdQ = c 2 / tcV0 • hdr (61) 

io Jfi Jr 

recall the definition of the velocity potential to transform the right-hand side, 


/ ivcpdtt + (? / VwV^dQ = c 2 / wu-ndr. (62) 

Jn Jn Jr 

Choosing Galerkin test functions, the derivatives of the trial functions with respect to the 
unknowns, the two volume integrals become 


[M,\ = / {i/} T {i/}dfi 

Jn 


[k,\ = / {vh} t {vh} an 

Jn 


where {H} is the vector of nodal interpolation functions. In the usual Finite Element (FE) 
matrix-vector form we get 


[Mf]{(j)} + c 2 [Kf]{(j)} = c 2 J wv ■ hdr. 
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Equations (62) and (65) hold for each element of a finite element domain, and, if 
continuity between elements is enforced, also hold globally. Thus the right-hand side is 
only defined on the boundary of the domain. At the fluid-structure interface the velocity 
compatibility provides a convenient input to this finite element problem. Specified Dirichlet 
or Neumann boundary conditions can also be applied with relative simplicity, but for this 
work non-reflecting boundary conditions are most appropriate yet are not easily applied. 

C. CELLULAR AUTOMATA MODEL OL THE WAVE EQUATION 

Cellular Automata (CA) are discrete, rule-based numerical methods that can model 
complex physical phenomena with relative simplicity. Generally, both space and time are 
treated discretely and the value of the quantity in question is limited to a finite set of val¬ 
ues. As the space-time domain proceeds or grows the seemingly simple model converges 
to the complex real-world behavior. The simplicity of the chosen rules and their imple¬ 
mentation lowers the computational cost while still achieving required accuracy. CA rules 
developed for modeling wave propagation are pre-cursors to the lattice Boltzmann method 
of modeling fluid flow. 

Following the work of Chopard [25, 26], Kwon and Hosoglu [28] modeled the wave 
equation in one- and two-dimensions with fixed boundaries using the following rules: 

<j>c{t + At) = - 4>c(t - At) + (f> E {t) (66) 

4>c(t + At) = (<pw{t) + 4> E (t ) + 4>s(t) + 0jv(t) - 2 4>c{t - At))/2. (67) 

The value of 0 at each interior grid point in the domain of interest (0c) is updated according 
to the values at its nearest Von Neumann neighbors (0 N , 05 , <f> E , <j) W ) as shown in Figure 
39. 

For convenience the set of points is divided into two sets, “black” and “white” (or 
“odd” and “even”), such that the neighbors of each white point are all black, and only one 
“color” is updated during each iteration. This model includes an evolution of the variable 0 
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Figure 39: Center node (black) and its Von Neumann neighbors (white) 

from being restricted to integer values (as in a traditional CA model) to being real valued. 
In its fully discrete form this CA rule was developed to model particle motion; with the 
relaxation that allows real-valued states, it also corresponds to the finite difference model 
of the wave equation on a uniform grid. 

1. One Dimension 

The classic illustration of D’Alembert’s solution to the one-dimensional wave equa¬ 
tion is the perturbation of a string subject to tension. For the moment, we shall apply fixed 
boundary conditions to the ends of the “string” and focus our attention to interior points 
well away from those ends. Further discussion of appropriate boundary conditions will 
follow. Consider a string subject to a Gaussian perturbation at its midpoint as shown in 
Figure 40, 


Utt = C 2 U XX — OO < X < oo 0 < t < oo 

_ A 

u(x, 0) = f(x ) = e 2 


( 68 ) 
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Mt(a,0) = g(x) = 0. 


The D’Alembert solution to this problem is 


1 (x — ct) 2 ( x+ct) 2 

u{x,t) —-[e 2 -+ e 2 


( 69 ) 


which agrees nicely with the CA solution computed using the rule found in Equation (66) 
as shown in Figure 41. 



Figure 40: Initial perturbation of an infinite string 
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Figure 41: CA solution vs. D’Alembert’s solution to the one-dimensional wave equation 
in an infinite string 

2. Boundary Conditions 

Fixed boundary conditions are easy to implement, but are of limited utility in mod¬ 
eling a potentially infinite fluid domain. One possible approach is to model a much larger 
domain than that of interest so that the area of interest is well away from the boundary and, 
as such, solutions within it are unpolluted by whatever boundary condition is imposed. This 
requires many computations that will be ignored-a seeming waste. Another approach is to 
generate a model of a non-reflecting boundary such that the wave in question is unaffected 
by its proximity. This is an active area of research in the finite element arena. In the CA 
arena, Chopard and Droz [25] suggested implementing boundary conditions by generat¬ 
ing virtual cells adjacent to the boundary cells as shown in Figure 42. Fixed or specified 
boundary conditions do not require a virtual neighbor, but are shown for completeness. 
Non-reflecting (or zero-gradient or adiabatic) boundaries model the behavior in the heart 
of the domain of interest, well away from the influence of any boundary. The free boundary 
condition corresponds to that of a constant gradient. In practice, this can be implemented 
either by generating the virtual neighbor cells or by applying the resulting rule to the actual 
cells located on the boundary in question. 
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fixed 


Figure 42: Virtual cell values for various boundary conditions in one dimension. After [25] 

Figure 43 shows the application of several different boundary conditions to the 
positive side of the spatial domain of our one-dimensional wave equation. For reference 
the f(x — ct) portion of the D’Alembert solution is shown as an exact solution. The non¬ 
reflecting boundary condition corresponds quite well with the analytic. All four waves peak 
initially in unison; the reflecting wave (black) returns with equal amplitude while the free 
wave (red) returns with an inverse amplitude. The fluid domain to be modeled will have 
non-reflecting boundary conditions imposed at the arbitrary edge of the domain and a free 
boundary used to represent the air-water interface. 

3. Discretization and Model Fidelity 

The various waves modeled above by CA rules display a coarseness that results 
from updating the value of each point in space at alternating time steps. The initial per¬ 
turbation displayed thus far has been a Gaussian wave of medium width that proved to 
be smooth enough to demonstrate the desired characteristics. Attempts to model a point 
source along the lines of f(pc) = 0, x ^ 0, f(x) — 1, x — 0 were unsuccessful, begging the 
question of how smooth a function or discretization are necessary to use CA to model the 
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traveling wave vs time at a common spatial point 



Figure 43: Application of various boundary conditions to CA calculation of one¬ 
dimensional wave propagation 


wave equation. In order to determine a sufficiently fine discretization relative to the sharp- 

/ -A 

ness of the initial perturbation, f(x) = e ^, a series of progressively narrower Gaussians 
as shown in Figure 44. For each initial condition, the CA solution to the wave equation was 
calculated for the same set of discretizations (progressively smaller dx) and error norms rel¬ 
ative to the D’Alembert solution were calculated at the same arbitrary time. Non-reflecting 
boundary conditions were applied in all cases. The resultant convergence as a function of 
dx is shown in Figure 45. An error norm of 1% was chosen as the comparison point. The 
largest dx value required to achieve that level of convergence was then plotted against the 
width at half maximum for each initial condition as shown in Figure 46. A linear estimate 
of that data is that dx c = |cr. Applying this estimate to the specific f(x) used, 4.5 elements 
or nodes are required to represent the descent from peak value of / to 1% of that peak. 
Therefore, a discretization that uses eleven or more nodes to represent both sides of a peak 
should be sufficient. 
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Figure 44: Initial perturbations of varying widths 



Figure 45: CA solution to Id wave equation convergence as a function of dx for various 
initial conditions 
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Figure 46: Critical discretization vs. initial perturbation width 
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4. Convergence 

To examine convergence of the present CA rule for wave propagation with respect 
to mesh or lattice spacing, comparison to the steady-state plane wave as presented by Junger 
and Feit [45] was used. Consider a semi-infinite fluid-filled space with a given uniform 
vibration, w{t) = We~ lult at the z = 0 boundary and a non-reflecting boundary as z 
oo. The steady state pressure in the wave guide is p(z,t ) = f)cWe t,kz ^ tuji> where k = 
The given function was applied to the z — 0 nodes in a CA domain with the initial 
values everywhere else uniformly zero, the CA rule was applied for a number of iterations 
corresponding to over three periods of the steady-state solution, and point-by-point error 
was calculated relative to the analytic solution over the range z G [0, and plotted in 
Figure 47. The rate of convergence is linear, as expected for a first-order method. 
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Figure 47: Convergence of CA wave equation rule to analytic solution as a function of dx 
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5. 


Three Dimensions 


It follows from Equations (66) and (67) that the three-dimensional wave equation 
can be modeled as 


<f>c(t + At) 



wit) + (f>E(t) + <t>s(t) + (pN it) + — 3 <j>c(t — At)] (70) 


with the subscripts on the new terms standing for “front” and “back,” respectively. To test 
this supposition, consider a “point” source located in the center of a domain of interest. 
As was discussed above, CA is not expected to faithfully model a true point source, so 
the source under consideration is a smooth radial function that has a maximum value of 
1 at the origin of the domain and is zero valued outside a radius of five nodes from the 
origin. The domain is comprised of 73 equi-spaced nodes in each direction; non-reflecting 
boundary conditions were applied on all six sides of the domain. Three points in space will 
be examined: the origin, an arbitrary point inside the initial perturbation (r < 5 dx), and 
an arbitrary point outside the initial perturbation (r > 5 dx). Referring to the notation of 
Equations (55) - (58), 


<p(x,y,z) 


1 — - if r < a, 
0 if r > a 


(71) 


fp(x,y,z) = 0 


(72) 


where r is the cartesian radius, \J x 2 + y 2 + z 2 . The analytic solution of the integrals in 
Equation (58) is not easily calculated, but application of a composite Simpson’s Rule over 
intervals of will yield a suitable comparison. 

Figures 48 - 50 show the analytical solutions at the respective points plotted as a 
function of time directly over their CA counterparts plotted versus the number of iterations 
through which the rule has been applied. These plots illustrate a key challenge to CA as 
identified by Hosoglu [46]: matching the discrete iterations of a cellular automaton to the 
continuous time domain, or calculating an appropriate time scaling factor (TSF). In one 
dimension, dt = — works well, but for three dimensions, as Figure 51 shows, there is a 
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phase difference resulting from a time scale mismatch. Consider a true point source located 
at the origin of an otherwise zero-valued CA domain. Following the current CA rule, the 
earliest a node at (dx, dy, dz ) can reach a value other than zero is after the third iteration; 
thus 3 dt = J dx 1 + dy 2 + dz 2 , or, for an equi-axed mesh, dt = c ^^ n - The points of 
interest for this exercise were chosen arbitrarily, but in such a way that both the inside and 
outside points are displaced from the origin in all three directions. Comparisons between 
the analytical and CA solutions plotted with this time equivalency are shown in Figures 
52-54. 


analytic solution at (0,0,0) 

1 ---.--- 


0 12 time 3 4 5 



Figure 48: Three-dimensional wave model at domain origin: time vs. iterations 
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Figure 49: Three-dimensional wave model at a point inside the initial perturbation: time 
vs. iterations 
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Figure 50: Three-dimensional wave model at a point outside the initial perturbation: time 
vs. iterations 
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Figure 51: Three-dimensional wave model at a point inside the initial perturbation: analytic 
solution vs time, CA solution vs — 

’ C 


origin 



Figure 52: Three-dimensional wave model at a domain origin 
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Figure 53: Three-dimensional wave model at a point inside the initial perturbation 
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Figure 54: Three-dimensional wave model at a point outside the initial perturbation 
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D. COUPLING OF FINITE ELEMENT AND CELLULAR AUTOMATA MOD¬ 
ELS OF THE WAVE EQUATION 

The CA model of the wave equation for velocity potential is simple to implement 
and can easily be adapted to a variety of non-trivial boundary conditions, but converting 
velocity potential back into useful quantities like pressure and velocity as time- and spatial- 
derivatives is hampered by the alternating update nature of the model. The finite element 
model, on the other hand updates the value of every point every time-step, but can be com¬ 
putationally expensive and non-trivial boundary conditions can be difficult to implement. 
A combination of the two methods would seem to resolve the short-comings of each and 
enable a faithful model of the fluid portion of our fluid-structure interaction. 

The general idea is to have several layers of finite elements in contact with the 
structure and then to have that fluid volume surrounded with a CA domain upon which the 
non-reflecting and free boundary conditions can be imposed as shown in Figure 55. The 
two fluid domains will overlap such that the outer layer of finite element nodes will be pro¬ 
cessed as interior CA nodes whose CA-calculated values become FE-specified boundary 
values. The next set of FE nodes inside the domain are calculated by the FE machinery and 
then passed to the CA domain to serve as neighbors for application of the CA rule to the 
outer set. These node sets are illustrated in Figure 56. 

The pre-validation of this scheme was to establish comparable domains of each 
model, impose a specified velocity potential field on one face of the domains (the top), 
and specify a fixed, zero-valued boundary condition on the other five faces. The specified 
input is a radially scaled sinusoid-it achieves its maximum value at the center of the face 
over which it is applied and is zero beyond a radius of one-half the width of the region. 
The resulting value of the velocity potential at the respective domain centers compares 
favorably as shown in Figure 57. 

The first full validation of the coupling scheme was to model a joint domain with 
an imposed velocity potential on the top of the FE portion of the domain which in turn 
rests atop the CA portion. The four sides of both domains have fixed zero-valued boundary 
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finite element domain within cellular automata domain 




y axis 


x axis 


Figure 55: Finite Element fluid domain surrounded by Cellular Automata fluid domain 

conditions and the bottom of the CA portion is non-reflecting. The value of the velocity po¬ 
tential on either side of the interface is shown in Figure 58. Next, the specified function on 
the top of the FE domain (the same radially scaled sinusoid) was treated as a velocity rather 
than velocity potential. Again, the values of </> are compared near the interface between the 
two models and shown in Figure 59. Next, the domain and interfaces are expanded such 
that the CA domain surrounds the FE domain on five sides and the non-reflecting boundary 
condition is applied to all six sides of the outer domain with the exception of the top of the 
FE domain-velocity is specified there. Those results are shown in Figure 60. 

The time-integration of Equation (65) was performed using a Newmark-/? algorithm 
with a zero valued [C] matrix.[47] The cases that involved fixed, zero boundary conditions 
(Figures 57 - 59) display noticeable high-frequency noise that is suspected to be caused 
by the sudden change in </> value at the boundaries. Despite a smooth input function, the 
beginning of a similar phenomenon is noticeable in the case with completely non-reflecting 
boundary conditions as well (Figure 60). The particular Newmark-/? scheme used (7 = |, 
(3 = |) is unconditionally stable and was chosen to ameliorate any potential difficulty 
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O FE calculation passed to CA 
+ CA calculation passed to FE 



Figure 56: Node sets for exchange of data between fluid domains 


arising from dt being a fixed function of dx, however a transition to an o-mcthod (or 
Hilber-Hughes-Taylor (HHT)) may be necessary to dampen this noise. See Appendix A 
for details on these methods. 
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velocity potential at center of domain 



Figure 57: Comparison of velocity potential propagation between finite element and cellu¬ 
lar automata models with common Dirichlet boundary conditions 


velocity potential near (0,0,z), <|>(x,y,^ o ,t) specified 



Figure 58: Velocity potential propagation between finite element and cellular automata 
domains, velocity potential (©(top)) specified 


71 









velocity potential near (0,0,z.), v(x,y,^ Qp ,t) specified 



Figure 59: Velocity potential propagation between finite element and cellular automata 
domains, velocity (v(top)) specified 


velocity potential near (0,0,z.), v(x,y,^ op ,t) specified 



Figure 60: Velocity potential propagation between finite element and cellular automata 
domains, FE inside CA, velocity (u(top)) specified, non-reflecting boundary conditions 
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1. Comparison with Homogeneous Fluid Domain 

Figure 61 demonstrates close agreement between the response of a composite do¬ 
main and a homogeneous fluid domain subject to the same input and boundary conditions. 


velocity potential at mid-domain 



Figure 61: Comparison of velocity potential at mid-domain resulting from specified value 
on one face: FE + CA domain vs. homogeneous CA domain 


73 






THIS PAGE INTENTIONALLY LEFT BLANK 


74 



V. RESULTS 


A. ACOUSTIC FIELD FLUID-STRUCTURE INTERACTION 

This chapter will combine the structural and fluid models discussed so far and ex¬ 
amine their interactions. Finally, comparison with experimental results will be presented. 

1. Information Exchange 

As previously discussed, the fluid model used in this work is that of the velocity 
potential. At each time step the transverse velocity of the structure is transmitted into the 
fluid where it is converted to velocity potential and propagated through the fluid domain ac¬ 
cording to the wave equation. The resulting fluid pressure, a function of the time derivative 
of the velocity potential is then applied to the structure via its load vector for calculation 
of displacement and velocity during the next time step. Unconditionally stable time inte¬ 
grators across the various domains allows the selection of time step size based on the CA 
time-scaling factor. 

For simplicity the meshes of the structure, the FE fluid, and the CA fluid are mutu¬ 
ally conforming. 

2. Homogeneous Isotropic Single-Layer Structure 

The algorithm described above was validated by modeling a clamped foot square 
plate one-quarter inch thick homogeneously comprised of an isotropic material. Both CG 
and DG models were generated. Fresh water material properties were used for the fluid 
model. Figure 62 shows the transverse displacement of the center of the plate subject to 
a constant concentrated force applied at its center. The dry case is the resulting oscil¬ 
lation about its predicted static deflection; the wet case shows that both magnitude and 
frequency of the oscillation have been altered as a result of the FSI. Kwon [48] discussed 
this phenomena and observed the effects of the elastic modulus and density of the struc¬ 
ture. Specifically, he noted that a structure with a density close to that of the fluid would 
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be more affected by the interaction. This observation is borne out in Figures 63 and 64 for 
which the respective material properties were changed from the baseline case. The baseline 
model has a structural density 2.7 times that of the fluid, a frequency ratio (dry/wet) of 1.84 
and an amplitude ratio (first peak - first trough, dry/wet) of 2.49. The double modulus case 
(Figure 63) shows a slightly higher frequency oscillation about a lesser static deflection in 
the dry case; the dry/wet frequency ratio is 1.86 and amplitude ratio is 2.15: a lesser relative 
change. The double density case (Figure 64) shows the expected lower frequency oscilla¬ 
tion in the dry case, but also a dry/wet frequency ratio of 1.49. The dry/wet amplitude ratio 
for the double density case is 2.00. 



Figure 62: Displacement of clamped plate with and without fluid-structure interaction 
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Figure 63: Displacement of clamped plate of double modulus with and without fluid 
structure interaction 



Figure 64: Displacement of clamped plate of double density with and without fluid 
structure interaction 
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3. Two-layer Plates 

The clamped plates modeled in this section are each 0.3048m x 0.3048m x 3.5mm 
and comprised of two thickness layers and sixteen elements in each planar direction. The 
material is an isotropic approximation of E-glass. Initial conditions were zero displacement 
and velocity; a constant concentrated force of 1000N applied at the center of the plate at 
the first time step. The “damaged” plates had a four element by four element debonding 
patch inserted between the two layers at the center of the plate. This debond was of the 
partial disconnection method described in Chapter III. Figures 65 - 67 display the time 
histories of displacement, velocity, and normal strain in the plane on the bottom of each 
plate. The stress profiles in Chapter III showed that maximum values were at the edges of 
the damage zone; similar phenomena are expected in strain values, but as Figure 68 shows, 
the strains at the center of the plate for both dry and wet cases is of greater magnitude than 
those at the +y edge of the damage zone (the maximum numerically of the four edges). 
This further suggests that an interface layer is necessary to properly model debonding in 
laminated composites and other layered structures. The density of E-glass is twice that of 
water, close enough to observe an added mass effect in the case of FSI. That is, a structure 
with density comparable to that of fluid reacts not only to the pressure effect of the fluid- 
dampening the amplitude of its vibrations, but also at a lower frequency as though it was a 
more massive structure. 
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damaged two layer E-glass plate displacement 



Figure 65: Displacement of damaged clamped two layer E-glass plate with and without 
fluid-structure interaction 


damaged two layer E-glass plate velocity 



Figure 66: Velocity of damaged clamped plate two layer E-glass plate with and without 
fluid-structure interaction 
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Figure 67: Strain at center of clamped two layer E-glass plate with and without fluid- 
structure interaction 



Figure 68: Strain at center and edge of damage zone of clamped two layer E-glass plate 
with and without fluid-structure interaction 
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4. Three-layer Plates 

The two-layer model of the previous sub-section did a poor job of reflecting debond¬ 
ing damage within a laminated plate. A three-layer model with a thin interface layer with 
properties approximating a common adhesive inserted between two layers of E-glass was 
subjected to the same loading and boundary conditions (zero initial displacement and ve¬ 
locity, 1000 N concentrated force at center, clamped edges); responses were calculated for 
five-hundred time steps. Figures 69 and 70 show that the displacement and velocity re¬ 
sponses of the center of the plate do not reflect the presence or absence of a debonding 
zone, but do demonstrate FSI mass effects in a fashion similar to that of the two-layer 
model. Figure 71 shows that the strain calculated in the E-glass elements reflects presence 
or absence of damage only mildly. Figure 72, on the other hand, shows clearly that the 
interface layer is profoundly affected by the presence of a damage zone. To be clear, the 
strain values at the centers of the interface layers of the damaged plates are not identically 
zero, but they are four orders of magnitude lower than their undamaged counterparts. 


three layer E-glass plate displacement 



Figure 69: Displacement of clamped three layer E-glass plate with and without fluid- 
structure interaction and with and without damage 
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three layer E-glass plate velocity 



Figure 70: Velocity of clamped three layer E-glass plate with and without fluid-structure 
interaction and with and without damage 



Figure 71: Strain of clamped three layer E-glass plate with and without fluid-structure 
interaction and with and without damage at center of lower E-glass layer 
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Figure 72: Strain of clamped three layer E-glass plate with and without fluid-structure 
interaction and with and without damage at center of interface layer 
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Figure 73: Strains in dry clamped three layer E-glass plate with and without damage 
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Figure 74: Strains in wet clamped three layer E-glass plate with and without damage 
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B. EXPERIMENTAL VALIDATION 


Recent experimental work by Kwon and his students [49], [50], [51], [29] examined 
the response of composite plates to low velocity impact with and without fluid-structure in¬ 
teraction^). In general, they found that for a given impact weight dropped from the same 
height, structures with lower density relative to the fluid in question experienced higher 
resultant forces and consequently greater damage than the same material in dry conditions. 
Additionally, they noted that the initial observable damage mode was delamination occur¬ 
ring on the face of the plate opposite the impact site. 


Carbon fiber layers 



Figure 75: Schematic of VARTM for plate manufacture. From [49] 


In his work, Conner [29] used a Vacuum Assisted Resin Transfer Molding technique, 
shown in Figure 75, to construct a series of 12in by 12in composite plates comprised of 
sixteen layers of E-glass (approximately 3.5mm thick in toto) and subjected them to low- 
velocity impact forces that resulted from dropping a 10.8kg weight from various heights to 
the center of the plate(s) using the assembly shown in Figures 76 and 77. The plates were 
instrumented with strain rosettes at four set positions and oriented such that one channel 
returned e x directly; e y was calculated as a function of all three channels and the included 
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Figure 76: Drop Weight Rig used in Impact Testing. From [48] 

angles as can be found in many solid mechanics textbooks including [52]. Data was sam¬ 
pled at a frequency of 10,000 Hz (dt = 10 -4 s). 

Numerical comparison with this experimental data was conducted using a DG struc¬ 
tural model comprised of a single layer of plate elements with a discretization of twelve 
elements in each planar direction. The overall structure has a length to thickness ratio of 
87:1 and each element has a length to thickness ratio of 7.3:1. In this model the material 
properties used are those of E-glass, but treated as an isotropic material-the Young’s mod¬ 
ulus along its fiber direction was taken in all three directions. The mass matrix is lumped 
and therefore, diagonal. The dry response was calculated using the a-method time inte¬ 
gration of the equation of motion for the plate. The wet, or FSI, response was calculated 
according to the acoustic field FSI described in the last chapter with time step size for the 
entire model equal to the TSF for the CA portion of the fluid model. 

The force inputs to the numerical plates were smoothed versions of the experimen¬ 
tally measured force data in the time region of interest-the main impact. The smoothing 
was conducted by sampling the raw data every five time steps, generating an interpolation 
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Figure 77: Drop Weight Rig as used for Impact Testing with FSI. From [48] 

function (using MATLAB’s interpl function with the spline option) and evaluating 
the interpolation function at every required time step for the model. For the dry cases, the 
same time step size as the experimental was used. Figure 78 shows the raw experimental 
force data and the smoothed version as generated for the dry case; Figure 82 shows the 
force data used for the FSI case. 

The calculated time history of the displacement field was used to calculate a time 
history of the strain vector at nodal points throughout the domain of the plate. Those nodes 
closest to the positions of the strain gages in the experimental work were examined relative 
to the recorded data and are plotted in Figures 79-81 for the dry plate and Figures 83 - 85 
for the FSI case. All show good qualitative agreement between experimental and numerical 
data. Differences can be attributed to the smoothing of the input force, ignoring impact 
effects, the isotropic treatment of an orthotropic material, the approximation of structural 
thickness, approximate positioning of strain gages, and mis-alignment of strain gages. The 
data correlating with gage 1 appears to show better overall agreement than the other two, 
most likely because that gage was located approximately equi-distant from both the point 
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of impact and the clamped boundaries of the plate. The other gages were closer to, and 
therefore, more exposed to the effects of the physical boundaries of the plate. 


force input 



Figure 78: Raw and smoothed experimental force data for dry plate 
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Figure 79: Measured versus calculated strain, dry plate, gage 1 
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Figure 80: Measured versus calculated strain, dry plate, gage 2 
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Figure 81: Measured versus calculated strain, dry plate, gage 3 


force input 



Figure 82: Raw and smoothed experimental force data for wet plate 
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Figure 83: Measured versus calculated strain, wet plate, gage 1 
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Figure 84: Measured versus calculated strain, wet plate, gage 2 
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Figure 85: Measured versus calculated strain, wet plate, gage 3 
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VI. CONCLUSION 


A. SUMMARY OF FINDINGS 

The goal of this work was to develop computational techniques to accurately model 
and simulate dynamic responses and failures of composite structures in an acoustic field. 
After implementing a nodal three-dimensional element to verify basic computational method¬ 
ology, a displacement-only plate finite element was formulated and implemented using 
Discontinuous Galerkin (DG) methodology. Such a displacement-only element allows con¬ 
struction of multi-layered structures like sandwich plates and other laminated composites 
in a manner similar to full three-dimensional solid finite elements. Results generated from 
this formulation compare favorably with theoretical predictions as well as existing CG nu¬ 
merical models for both static and dynamic responses for both simple and multi-layered 
plate structures. 

Application of the new element to the analysis of failure initiation and propagation 
in sandwich composite structures shows great promise. Static qualitative stress profiles are 
similar to those found using CG techniques, but the elemental rather than nodal connectivity 
used in DG formulations suggests a simple means of modeling debonding between material 
layers by disconnecting their respective elements in the global stiffness matrix. Complete 
disconnection of neighboring elements in an imposed debonding zone was shown to be 
incorrect because it allowed the core layer to deflect not only through the disconnected 
resin layer but also through the still-present skin layer. Partial disconnection-removing 
connectivity between opposing pairs of dofs in the planar directions but retaining weak 
connectivity for transverse pairs of dofs resulted in a stress profile that makes good quali¬ 
tative sense. Maximum stress values in both skin and core layers decreased in magnitude 
and moved from the center of the plate to the edges of the debonding zone, behaving like a 
stress concentration. The simplicity of this partial disconnection method can be a tremen- 
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dous computational savings when modeling the progression of damage without need for 
re-meshing or recalculation of the global stiffness matrix. 

Examination of FSI for the impact problem does not require a full fluid-flow model- 
the propagation of velocity potential according to the wave equation in the acoustic domain 
is sufficient for these purposes; as such, an extension of Cellular Automata (CA) from two 
to three-dimensions in modeling the acoustic field was demonstrated and validated. CA 
was chosen for this application not only because of the simplicity of its update rule, but 
also because of its flexibility in the implementation of non-reflecting boundary conditions. 
The alternating update nature of CA makes calculation of both spatial and time derivatives- 
required to convert the velocity of the structure into velocity potential in the fluid domain 
and to convert velocity potential into a pressure field-difficult. Insertion of a small finite el¬ 
ement interface zone between the structure and the CA fluid domain resolved this difficulty 
for a relatively low computational cost. The combination of a small FE acoustic domain 
with an enveloping CA domain proved to be an efficient way to implement non-reflecting 
boundary conditions. 

Finally, the combined model of a DG structure interacting with a FE+CA fluid do¬ 
main was shown to have good agreement between calculated and experimentally measured 
strain values for plates subject to low-velocity impact in the pre-damage regime. In partic¬ 
ular, the added mass effect on structures with low density relative to the fluid medium was 
apparent in both simulation and experimental comparisons. 

The methods developed and examined in this study: the displacement-only DG 
plate finite element, the partial disconnection failure model, and the hybrid FE+CA acoustic 
field model, show great promise for flexible and accurate modeling of debonding initiation 
and propagation in sandwich and laminate composite structures subject to FSI. 

B. FUTURE WORK 

This work should be viewed as a starting point for further investigation into the 
utility of DG methods in composite failure modeling. 
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First and foremost, the current formulation must be re-implemented to be a true 
element-wise computation in order to reap the benefits of DG not just pay the costs of solv¬ 
ing for a greater number of dof. This type of update will enable more efficient and flexible 
modeling of larger problems, including more refined meshes, more complex geometry, and 
approaches that address all levels of multi-scale analysis of composite materials. Such a 
re-implementation should also enable the coupling of the current DG formulation with CG 
codes. This should be readily achievable as the current formulation is derived from one 
with such coupling as a specific goal. This sort of coupling can be used as an alternative to 
refining the mesh in the areas of interest like existing or expected damage by replacing the 
refined mesh with DG elements in order to better examine the physical phenomena. 

A more computationally efficient implementation should also include and enable 
progressive failure modeling both through inclusion of traction-separation type models for 
the post-damage regime, but also through propagation of damage beyond its initiation site. 
Addition of a full impact-impulse model for force input should enable even more faithful 
modeling and closer comparison with experimental results. 

The closer a computational model approaches observed physical phenomena, the 
more useful and trustworthy its results in evaluating more complex geometries and operat¬ 
ing environments. 
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APPENDIX A. TIME INTEGRATION ALGORITHMS 


This appendix contains explication of the algorithms used to solve the matrix-vector 
equation of motion in this work. Implicit methods are favored for their stability independent 
of size of time step-a concern when matching various domains. All methods are trying to 
solve 

[M](u} + [ C\{u} + [AIM = (/} (73) 

for {u}. Solutions for the two time derivatives are used as needed to update {»}; in the 
finite element fluid domain {-u} is also used for a pressure calculation. 

A. NEWMARK-/3 METHOD 

The below algorithm, taken from [47] was initially implemented to serve as time 
integrator of Equation (73) by using the weighted averages 

{u } n +1 = {u}n + [(1 - l){u} n + l{u}n+l] ' dt (74) 

dt 2 

{u}n +1 = {u}n + {u}n ■ dt + [(1 - 2 fd){u} n + 2(3{u} n+1 \ ■ — (75) 

substituting and rearranging terms results in 

[M + ^dtC + (3dt 2 K}{u} n+1 = {f} n+1 

dt 2 

- [(1 - 7 )dtC + (1 - 2 j3)—K]{u} n - [C + dtK]{u} n - K{u} n (76) 

which can be solved for u and then u and, in turn, u at each time step. The method is 
unconditionally stable for 2f3 > 7 > |. Parameter choices are 7 = \ and /3 = \ correspond 
to Newmark’s Constant Average Acceleration Method. 

Dirichlet and Neumann boundary conditions can be imposed nodally by solving 
for the current acceleration on the boundary through Equations (74) and (75), zeroing the 
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corresponding rows of the compound left-hand-side matrix, setting the diagonal elements 
of those rows to 1 and substituting the boundary accelerations into the right-hand-side 
vector. 

B. q-METHOD 

Hilber, Hughes, and Taylor improved upon Newmark-/? with their introduction of 
the a-method [53]-[54]. This method is designed to dissipate high frequency noise without 
degrading the order of solution accuracy. The update rules for { «} and { «} are the same as 
in (75) and (74), but now the equation of motion is also a weighted average: 

[M]{u } n+ 1 + (1 + a)[C]{ii } n+ 1 - a[C]{u} n 

+ (1 + a)[K}{u } n+ 1 - a[K]{u} n = {f{t n+i+a )} (77) 


which is rearranged to 

[M + (1 + a)dt{jC + /3dtK)\{u} n+1 = (1 + a)f n+1 - af n 

-(l + a )*[(l- 7 )C + |(l-2/?)/v]{«}„ 

- [C + dt( 1 + a)K]{u} n - K{u} n (78) 

and then solved for («} n+1 which is then used to update {«} n+1 , and {u} n+ i. Dirichlet 
boundary conditions are applied nodally via rearrangement of Equation (75) to solve for 
prescribed values of the right hand side of Equation (78). This method is unconditionally 
stable when a £ [—|, 0], 7 = (1 — 2a)/2, and — (1 — cc) 2 /4. When a = 0 this method 
reduces to Newmark’s Constant Average Acceleration Method. 

C. TIME DISCONTINUOUS GALERKIN METHOD 

Another time integrator considered but not fully implemented in this work is the 
Time Discontinuous Galerkin (TDG) method presented by Chien, Yang, and Tang [55]. 
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They present time as yet another domain that can be discretized by finite elements, in par¬ 
ticular as discontinuous finite elements, as shown in Figure 86. For the undamped equation 
of motion ([C]=0 in Equation (73)), they generate the following matrix equation to solve 
for the displacement and velocity at each end of a particular time interval (or element): 
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where K and M are the usual stiffness and mass matrices and 


F\ = 



4>iit)Fdt 


( 80 ) 


F 2 = [ <j> 2 (t)Fdt (81) 

Jin 

where F is the usual applied load vector and 0 \it) and <f> 2 (t) are the time shape functions 
shown in Figure 86. Simplifying and recasting some terms results in 
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which can be solved for the velocity terms from which displacements may be calculated 
directly. They also show that this method is also unconditionally stable, making it an 
intriguing avenue for future work in structural dynamics. 



Figure 86: Temporal elements for TDG method. From [55] 
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APPENDIX B. IMPLEMENTATION DETAILS 


The discussion of the fluid and structural and fluid models in the main body of this 
thesis is (moderately) general and symbolic. All implementation was in MATLAB, and run 
on a MacBook Air (mid-2009), MacBook Pro (mid-2011), or on the Hamming cluster. 

A. MESHING 

The various domains were meshed for human convenience rather than any matrix 
bandwidth considerations. In general, node and element numbering proceeds from (x mm , 
Umin, Zrm.n) along the x-axis, then increment up the y-axis, and then up the z-axis. While the 
code does execute full and proper Jacobian calculations for arbitrarily oriented elements, in 
general the x-axis corresponds to the canonical r-axis; the y-axis corresponds to the canon¬ 
ical s-axis, and the z-axis, which also generally corresponds to thickness, corresponds to 
the canonical t-axis. 

The various meshes employed are conforming for human convenience. The fluid 
domain is entirely equi-axed. The FE portion is constructed using the coordinates of the 
bottom of the plate elements as a foundation and extending down a specified number of 
layers in the —z direction. The CA portion is constructed by specifying a factor by which 
plate length is multiplied-the cube of this value is the fluid volume. The CA nodes that are 
wholly inside the FE fluid domain and not needed are simply removed from the index sets 
and ignored. 

B. NUMERICAL INTEGRATION 

Exact integration is performed using Gauss-Lobatto quadrature. Under-integration 
of the shear terms of plate elements is performed using Gauss-Legendre quadrature. Be¬ 
cause the integration points for Gauss-Lobatto quadrature are still nodal (interpolation) 
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points, terms in K s that are products of different cardinal interpolation functions will be 
uniformly zero, resulting in a too sparse shear stiffness matrix. 

A numerical experiment comparing the static deflections calculated using the two 
different quadrature rules for calculating an under-integrated K s with both theoretical static 
deflection and that resulting from an exactly integrated K s was conducted. For a clamped 
plate of dimensions 0.3048m x 0.3048m x 0.00635m (12in x 12in x l/4in), elastic modulus 
70 GPa, Poisson’s ratio 0.3 subjected to a concentrated load of 1000N, predicted static 
deflection of the center of the plate is 0.31697mm [33]. Calculated deflection for Gauss- 
Lobatto quadrature was 0.0945mm; for Gauss-Legendre quadrature: 0.3086mm; for exact 
integration of K s : 0.0325mm. Clearly, Gauss-Lengendre quadrature is a better choice for 
under-integration of the shear stiffness matrix in this application. 

C. APPLICATION OF BOUNDARY CONDITIONS AND EXTERNAL LOADS 

Boundary conditions and external loads are applied nodally. For boundary condi¬ 
tions, rows of the mass and stiffness matrices corresponding to constrained dofs are zeroed; 
their diagonal elements are set to 1, and the corresponding element in the right-hand-side 
vector is set to the constrained value. External loads must first be converted to units of 
force and then distributed to the appropriate elements of the right-hand-side vector. For 
the pressure exerted by the fluid domain on the wet side of our notional plate, the nodal 
pressure vector is multiplied by a two-dimensional interpolation (mass) matrix to convert it 
to a force vector. For application to the DG structure, that resulting vector is decomposed 
to reflect the number of dof found at each geometric position; in this way, the force applied 
at a point ’’shared” by four separate discontinuous elements will be parsed among them 
equally. 
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D. MATLAB SPECIFICS 


1. Sparse Matrices 

The lumped mass and global stiffness matrices are sparse, the former is diagonal 
and the latter is block tri-diagonal. MATLAB does allow assembly of a full stiffness matrix 
followed by K=sparse (K), but the calculation and assembly designed from the outset to 
be sparse is a more efficient use of computing resources. To do so, the global indices of 
each element of the collection of elemental K matrices are stored in vectors i v and j v 
with the corresponding values stored in k v. After the elemental calculations are complete, 
K=sparse (i_v, j_v, k_v) returns the sparse global stiffness matrix. 

2. CA Implementation 

The CA portion of the fluid model is updated using sets of node indices. The domain 
is comprised of N nodes where N — N x x N y x N z . In this work N x — N y — N z =an 
odd number; this keeps the eight corners of the domain in the odd set and is convenient, 
not required. The CA coordinate array is used to match node numbers of corresponding 
geometric points between the two fluid domains; an N x 6 array named neighbor tracks 
the eponymous relations by node number with a 0 entry indicating the end of the domain 
in that direction, the velocity potential in CA domain is described by the two A-vectors, 
phi and phi old. A large portion of the CA set-up is the definition of index sets for these 
two vectors. The largest two are the odd and even interior points, oint and eint. The six 
faces and twelve edges of the rectangular fluid domain are likewise split into odd and even 
index sets. This infrastructure makes an iteration of the CA rule a simple matter of calling 
subroutines UpdateFace or UpdateEdge with arguments of phiold, neighbor, 
the index set of the area to be updated, and a flag indicating the boundary condition to be 
employed. 
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E. MATERIAL PROPERTIES 

The following tables specify the material properties used for various calculations in 
this study. 


Property 

Value 

G xz 

345 MPa 

r< 

^Tyz 

345 MPa 


Table 2: Material Properties of aluminum honeycomb. From [41], [42] 


Property 

Value 

E 

72.4 GPa 

V 

0.3 


Table 3: Material Properties of aluminum skins, From [41], [42] 


Property 

Value 

Ex 

17.24 GPa 

P 

2020 

ITT 5 

L'xy 

0.3 

G xy 

6.619 GPa 

Ey 

17.24 GPa 

E z 

7.929 GPa 


0.24 

Vyz 

0.24 

Gxz 

2.896 GPa 

r 

^Tyz 

2.896 GPa 


Table 4: Material Properties of E-glass, From [56] 
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Property 

Value 

E 

8.34 GPa 

P 

1180^1 

in ' 3 

V 

0.28 


Table 5: Material Properties of Epoxy Resin, From [57] 
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